rm(list=ls())
###Monte Carlo Code for ECM
library(lattice)
library(AER)
library(matrixStats)
library(stats)
#Open Text File for Results
#sink("c:/test")

##########################
# Lag Operator Funcction #
##########################

tslag<-function(y,d=1){
n<-length(y)
c(rep(NA,d),y)[1:n]
}

#
# Generate an autoregressive X 
# rho gives the degree of autocorrelation in x
# n is the samle size
#   

gen.x <- function(n, rho){
	e <- rnorm(n)
	x <- double(n)
	x[1] <- rnorm(1)
	for(i in 2:n){
		x[i] <- rho * x[i-1] + e[i]
	}
	return(x)
}

# 
# Create a function to generate Y as an ADL(1,1)
#
#   n is sample size
#   a and r are the short run effects of x and lag x
#   g is the dynamic effect of y
#   x is the weakly exogenous covariate
#

gen.y <- function(n,a,g,r,x){
	e <- rnorm(n)
	y <- double(n)
	y[1] <- rnorm(1)
	for(i in 2:n){
		y[i] <- g*y[i-1] + a*x[i] + r*x[i-1] + e[i]
	}
	return(y)
}


#
# Function to generate the data frame, estimate the model, and save 
# the model and  the model parameters. 
#
# The function takes as inputs:
#   Coefficients
#   Sample Size
#   And the Number of Montecarlo replications
#   
# The function returns a list of the results for each replication,
# and the experimental values (a,r,g,rho, and n).
#

a<- .25
r<- .5
g<- .5

experiment <- function(a=.25, r=.5, g=.5, rho=0, n=1000, nsim=1000){
  #  set.seed(j)
  res <- vector(mode="list", length=nsim)
  for(m in 1:nsim){
    
    # generate x
    x <- gen.x(n=n, rho=rho)
    # generate y
    y <- gen.y(n=n, a=a, g=g, r=r, x=x)
    # generate lags and differences in x
    x.diff <- x - tslag(x)
    y.diff <- y - tslag(y)
    lagx <- tslag(x)
    lagy <- tslag(y)
    
    #Create Data Frame
    simecm <- data.frame(y, x, x.diff, y.diff, lagx, lagy)
    # run ADL model
    adl <- summary(lm(y~ lagy + x + lagx,data =simecm))
    # run restricted ADL model
    adl.r <- summary(lm(y~ lagy,data =simecm))
    # F test for ADL
    aov.adl <- anova(lm(y~ lagy + x + lagx,data =simecm),lm(y~ lagy,data =simecm))
    adl.F <- aov.adl$F[2]
    f.crit <- qf(.95, df1=2, df2=n-4) 
    adl.F.infer <- ifelse((adl.F)>f.crit,1,0)   
   
    #LRM ADL
    
    b.adl <- adl$coef[,1]
    bt.adl <- t(b.adl)
    lrm.adl.est <- ((b.adl[3] + b.adl[4])/(b.adl[2] - 1))
    
    # run error correction model
    ecm <- summary(lm(y.diff ~ lagy + x.diff  + lagx, data = simecm))
    # save model and parameters 
    # Comparing the ADL and ECM Inferences
    ecm.t.lagx <- ecm$coef[,3][4]
    t.crit <- abs(qt(.025,n-4))
    ecm.t.infer <- ifelse(abs(ecm.t.lagx)>t.crit,1,0)
    ecm.adl.x.compare <- ifelse(adl.F.infer == ecm.t.infer,1,0)
    #IV Reg for LRM SE
    inst <- cbind(1,x, lagy, lagx)
    ivr <- summary(ivreg(y ~ y.diff + x + x.diff | x + inst, data = simecm))
    # Bias
    
    true<-c(0,(g-1),a, (a+r), ((a+r)/(g-1)))
    
    b.ec <- cbind(ecm$coef[,1])
    bt.ec <- t(b.ec)
    estses.ec <- cbind(ecm$coef[,2])
    estparms.ec <-matrix(c(b.ec,b.ec[4,]/b.ec[2,]))
    bias.ec <- matrix(estparms.ec - true)
    
    # MSE
    
    par <- c(0,(g-1),a, (a+r))
    mse.ec <- matrix((bt.ec-par)^2)

    # Results
    
    res[[m]] <- list(mod = ecm, mod.adl = adl,rmod.adl = adl.r,aov.adl=aov.adl, 
                     adl.F = adl.F, adl.F.infer=adl.F.infer,ecm.t.lagx=ecm.t.lagx,
                     ecm.t.infer=ecm.t.infer,ecm.adl.x.compare=ecm.adl.x.compare,
                     a=a, r=r, g=g, rho=rho, n=n, mod1 = ivr,estses.ec=estses.ec,
                     bias.ec=bias.ec,mse.ec=mse.ec,lrm.adl.est=lrm.adl.est)
  }
  class(res) <- "ecmsim"
  return(res)
}



#
# Run the experiments
#
#   (n) is a list of values for the different sample sizes you want to
#   calculate for each rho.
#
#   (rho) is a list of autoregressive parameters for X
#
#   (g) is the list of autoregressive parameters for y 
#
#   (a) the relationship between contemporaneous x and y
#
#   (r) the relationship between lag x and y
#
#   (eg) combines n and rho into a matrix that defines the dimensions of the
#   experiments and can be used as a reference. 
#
#   (allexp) holds all the results from all of the regressions for the
#   experiments. 
#

n <- c(25,50,75,100,250,500,750,1000)
rho <- c(.9)
g <- c(.5, .7, .9, .95, .98)
a <- c(.25)
r <- c(.5)
  
eg <- expand.grid(n, rho,g,a,r)
allexp <- vector(mode="list", length=nrow(eg))

system.time(

for(v in 1:nrow(eg)){
allexp[[v]] <- experiment(n=eg[v,1], rho=eg[v,2],g=eg[v,3],
                          a=eg[v,4],r=eg[v,5], nsim=1000)
}

)
#
# A function to summarize the results of the experiments
#
#   (par) is the vector containing the ture values of the parameters
#
#   (true) 

summary.ecmsim <- function(obj, par=c(0,(g-1),a, (a+r)), digits=3){
  
  
  ##
  ## ECM 
  ##
  
  ## GECM
  
  b <- sapply(obj, function(x)x$mod$coef[,1])
  bt <- t(b[1:4,])
	estparms <- cbind(t(b[1:4,]), b[4,]/-b[2,]) # 4 coefficients and LRM
	estses <- t(sapply(obj, function(x)x$mod$coef[,2])) # 4 standard errors for coefficients
	
  # GECM LRM
  
  lrm.ecm <- b[4,]/-b[2,]
  
  ## Bewley 
  
  b2 <- sapply(obj, function(x)x$mod1$coef[,1])
  bt2 <- t(b2[1:4,])
  estparmsivr <- cbind(t(b2[1:4,])) # 4 coefficients and LRM
  
  # Bewley LRM
  
  estlrm <- cbind(b2[3,]) # 4 coefficients and LRM
  
  # SE for LRM
  
  estsesivr <- t(sapply(obj, function(x)x$mod1$coef[,2]))
  lrmse <- cbind(estsesivr[,3])
  
  ##
  ## Reported Values
  ##
  
  #
  # t statistics
  #
  
  # LAG Y and LAG X from GECM
  
  t.na2 <- estparms[,2]/estses[,2] # t-statistic for error correction coefficients
  t.na4 <- estparms[,4]/estses[,4] # t-statistic for lagX
  
  # LRM Bewley
  
  t.LRM <- estparmsivr[,3]/estsesivr[,3]
  
  # Bias EC
  
  bias.sum <- sapply(obj, function(x)x$bias.ec)
  bias <- matrix(rowMeans(bias.sum))
  
  # MSE EC
  
  mse.sum <- sapply(obj, function(x)x$mse.ec)
  mse <- matrix(rowMeans(mse.sum))
  
  # Rejection Rates EC
  
	rejrateEC <- mean(ifelse(t.na2<(-1.645), 1, 0))
  rejrateLAGX <- mean(ifelse(abs(t.na4)>1.96,1,0))
  rejrateLRM <- mean(ifelse(abs(t.LRM)>1.96,1,0))

	cis <- estparms[,4] + cbind(-1.96*estses[,4], 1.96*estses[,4])
	covg <- mean(cis[,1] < par[4] & par[4] < cis[,2])
  estLRM <- estlrm
  LRMse <- lrmse
  
  ##
  ## ADL
  ##
  
  ## Parameters and Standard Errors
  
  b.adl <- sapply(obj, function(x)x$mod.adl$coef[,1]) # 4 coefficients
  estparms.adl <- cbind(t(b.adl[1:4,]))
  estses.adl <- t(sapply(obj, function(x)x$mod.adl$coef[,2])) # 4 standard errors for coefficients
  
  ## T Statistics
  
  t.LAGY.adl <- estparms.adl[,2]/estses.adl[,2] # t-statistic for lag y adl
  t.x.adl    <- estparms.adl[,3]/estses.adl[,3] # t-statistic for x adl
  t.LAGX.adl <- estparms.adl[,4]/estses.adl[,4] # t-statistic for lag x adl
  
  ## LRM ADL
  
  lrm.adl <- cbind((b.adl[3,] + b.adl[4,]) / (1 - b.adl[2,]))
  
  # Rejection rates for ADL
  
  rejrateLAGY.adl <- mean(ifelse(abs(t.LAGY.adl)>(1.96), 1, 0))
  rejrateX.adl <- mean(ifelse(abs(t.x.adl)>1.96,1,0))
  rejrateLAGX.adl <- mean(ifelse(abs(t.LAGX.adl)>1.96,1,0))
  
  ##
  ## Comparisons
  ##
  
  #
  # F test X and LagX in ADL vs. t test for LagX in ECM
  #
  
  # Note: A lot of this work is done in the experiment above.
  
  F.X.adl <- mean(sapply(obj,function(x)x$adl.F))
  
  f.x.adl <- sapply(obj,function(x)x$adl.F)
  samerate <- mean(sapply(obj, function(x)x$ecm.adl.x.compare))
  
  ft.bind <- cbind(t.na4,F.X.adl)
  adlVecm <- matrix(colMeans(ft.bind), ncol=1)
  
  
  # Compare LRMS
  
  LRMS <- cbind(lrm.adl,lrm.ecm,estlrm)
  colnames(LRMS) <- c("ADL","ECM","BEWLEY")
  
  # Organize output
  
  rownames(bias) <- c("Constant", "Error Correction Coefficient", 
                      "Contemporaneous Effect of X", "Lag X Effect", "LRM")
  rownames(mse) <- c("Constant", "Error Correction Coefficient", 
                     "Contemporaneous Effect of X", "Lag X Effect")
  colnames(mse) <- "MSE"
  colnames(bias) <- "Average Bias"
  cat("MSE:\n")
  print(round(mse, digits))
  cat("\nBias:\n")
  print(round(bias, digits))
  cat("\nAggregate MSE: ", round(mean(mse), digits), "\n")
  cat("Rejection Rate EC: ", round(rejrateEC, digits), "\n")
  cat("Rejection Rate LAGX: ", round(rejrateLAGX, digits), "\n")
  cat("Rejection Rate LRM: ", round(rejrateLRM, digits), "\n")
  cat("Coverage Rate: ", round(covg, digits), "\n")
  cat("LRM: ", round(mean(estLRM),digits), "\n")
  cat("LRMse: ", round(mean(LRMse),digits), "\n")
  
  
  rownames(adlVecm) <- c("t.lagX.ECM","f.XlagX.adl")
  colnames(adlVecm) <- "Test Statistics"
  
  invisible(list(bias=bias, agmse=mean(mse), mse=mse, rejrateEC=rejrateEC,
                 rejrateLAGX=rejrateLAGX,rejrateLRM=rejrateLRM, 
                 coverage=covg, tstat.ecc.ecm=t.na2,tstat.lagx.ecm = t.na4,
                 f.x.adl=f.x.adl, adlVecm=adlVecm,samerate=samerate,
                 LRMS=LRMS))
}


#
# Returns Elements of summary.ecmsim and t-stats
#

lapply(allexp,summary)


#
# Returns only the elements of sumary.ecmsim
#

for(i in 1:length(allexp)){
  summary(allexp[[i]])
}

#
# Look at individual elements of Allexp
#

# Look through 1-72 (9 rhos and 8 sample sizes)

s <- summary(allexp[[24]])
s








###############################################################
### Tables Showing Bias, Rejection, Coverage, and Inference ###
###############################################################

#
# Look at the Bias
#


bias <- round(sapply(allexp, function(x)summary(x)$bias),5)
bias.dat <- cbind(eg, t(bias))
colnames(bias.dat) <- c("T","RHO","G","A","R","CONSTANT","ECC","X","LAGX","LRM")
bias.dat

#
# Look at the MSE
#

mse <- round(sapply(allexp, function(x)summary(x)$mse),5)
mse.dat <- cbind(eg, t(mse))
colnames(mse.dat) <- c("T","RHO","G","A","R","CONSTANT","ECC","X","LAGX","LRM")
mse.dat

#
# Rejection Rates EC
#

rejectionEC.rate <- round(sapply(allexp, function(x)summary(x)$rejrateEC),5)
rejectionEC.dat <- cbind(eg,rejectionEC.rate)
rejectionEC.dat

#
# Rejection Rates LAG X
#

rejectionLAGX.rate <- round(sapply(allexp, function(x)summary(x)$rejrateLAGX),5)
rejectionLAGX.dat <- cbind(eg,rejectionLAGX.rate)
rejectionLAGX.dat

#
# Rejection Rates LRM
#

rejectionLRM.rate <- round(sapply(allexp, function(x)summary(x)$rejrateLRM),5)
rejectionLRM.dat <- cbind(eg,rejectionLRM.rate)
rejectionLRM.dat


#
# Coverage Rates
#

coverage.rate <- round(sapply(allexp, function(x)summary(x)$coverage),5)
coverage.dat <- cbind(eg,coverage.rate)
coverage.dat

#
# Matrix of Rates
#

output.mat <- cbind(eg,rejectionEC.rate,rejectionLAGX.rate,rejectionLRM.rate,coverage.rate)
colnames(output.mat) <- c("T","RHO","g","a","r","RR.ECM","RR.LAGX","RR.LRM","Cov.LAGX")
output.mat

#
# Comparing ECM and ADL
#

rate <- round(sapply(allexp, function(x)summary(x)$samerate),1)
comparison <- round(sapply(allexp, function(x)summary(x)$adlVecm),2) 
comparison.dat <- cbind(eg,t(comparison),rate)
colnames(comparison.dat) <-c("T","RHO","g","a","r","t.lagX.ecm","F.XlagX.adl","SameInferenceTandF?")
comparison.dat


sec2.1.dat <- cbind(eg[1:3],rate,bias.dat[,10],rejectionLRM.rate)
sec2.1.dat



#########################################
### Distributions of Test Statisitics ###
#########################################

#
# t-Statistics
#

# Pull t-Statistics out of Allexp

tstats.ecc <- sapply(allexp, function(x)summary(x)$tstat.ecc.ecm)
tstats.lagx <- sapply(allexp, function(x)summary(x)$tstat.lagx.ecm)
fstats.xlagx.adl <- sapply(allexp, function(x)summary(x)$f.x.adl)

# Calculate the Mean and Standard Deviations of the t-statistics
avg.tstats.ecc  <- colMeans(tstats.ecc)
sd.tstats.ecc  <- colSds(tstats.ecc)

avg.tstats.lagx <- colMeans(tstats.lagx)
sd.tstats.lagx <- colSds(tstats.lagx)

avg.fstats <- colMeans(fstats.xlagx.adl)
sd.fstats <- colSds(fstats.xlagx.adl)

# Calculate interesting quantiles 

t.stats.ecc.quantiles <- sapply(allexp, function(x)quantile(summary(x)$tstat.ecc.ecm, c(.025,.25,.5,.75,.975)))
t.stats.lagx.quantiles <- sapply(allexp, function(x)quantile(summary(x)$tstat.lagx.ecm, c(.025,.25,.5,.75,.975)))
f.stats.xlagx.quantiles <- sapply(allexp, function(x)quantile(summary(x)$f.x.adl, c(.025,.25,.5,.75,.975)))

# Bind values into a matrix

tstats.ecc.dat <- cbind(eg,avg.tstats.ecc,sd.tstats.ecc,t(t.stats.ecc.quantiles))
tstats.ecc.dat

tstats.lagx.dat <- cbind(eg,avg.tstats.lagx,sd.tstats.lagx,t(t.stats.lagx.quantiles))
tstats.lagx.dat

fstats.dat <- cbind(eg,avg.fstats,sd.fstats,t(f.stats.xlagx.quantiles))
fstats.dat

# Creates a dataset containing the tstats with corresponding n and rho

newdata.ecc <- data.frame(
  ststats = c(tstats.ecc), 
  n = rep(eg[,1], each=nrow(tstats.ecc)), 
  rho = rep(eg[,2], each=nrow(tstats.ecc))
)

newdata.lagx <- data.frame(
  ststats = c(tstats.lagx), 
  n = rep(eg[,1], each=nrow(tstats.lagx)), 
  rho = rep(eg[,2], each=nrow(tstats.lagx))
)

newdata.fstats <- data.frame(
  ststats = c(fstats.xlagx.adl), 
  n = rep(eg[,1], each=nrow(fstats.xlagx.adl)), 
  rho = rep(eg[,2], each=nrow(fstats.xlagx.adl))
)

#
# Plots t stats of the k experiments over rho
#


densityplot( ~ ststats | rho, groups=n, data=newdata.ecc)

#
# Plots the t stats for a given n over rho
#

densityplot( ~ ststats | n, groups=rho, data=newdata[which(newdata$n == 25), ])
densityplot( ~ ststats | n, groups=rho, data=newdata[which(newdata$n == 50), ])
densityplot( ~ ststats | n, groups=rho, data=newdata[which(newdata$n == 75), ])
densityplot( ~ ststats | n, groups=rho, data=newdata[which(newdata$n == 100), ])
densityplot( ~ ststats | n, groups=rho, data=newdata[which(newdata$n == 250), ])
densityplot( ~ ststats | n, groups=rho, data=newdata[which(newdata$n == 500), ])
densityplot( ~ ststats | n, groups=rho, data=newdata[which(newdata$n == 750), ])
densityplot( ~ ststats | n, groups=rho, data=newdata[which(newdata$n == 1000), ])

#############################################################################################

#############################################################################################
